Journal  of  Power  Sources  218  (2012)  192-203 


ELSEVIER 


Contents  lists  available  at  SciVerse  ScienceDirect 

Journal  of  Power  Sources 

journal  homepage:  www.elsevier.com/locate/jpowsour 


pri 

Sbb..«jtS 


Three-dimensional  modeling  of  polarization  characteristics  in  molten  carbonate 
fuel  cells  using  peroxide  and  superoxide  mechanisms 

M.Y.  Ramandi3’*,  P.  Bergb,  I.  Dincer3 

a Faculty  of  Engineering  and  Applied  Science,  University  of  Ontario  Institute  of  Technology,  Oshawa,  Canada  LIH  7I<4 
b  Faculty  of  Science,  University  of  Ontario  Institute  of  Technology,  Oshawa,  Canada  LIH  7I<4 


HIGHLIGHTS 


►  A  detailed  three-dimensional,  mathematical  model  of  an  MCFC  was  presented. 

►  Both  peroxide  and  superoxide  mechanisms  predicted  the  linear  region  of  the  polarization  curve  accurately. 

►  Both  mechanisms  showed  a  concave  tendency  for  the  cathode  reaction  rate  as  the  carbon-dioxide  mole  fraction  is  decreased  when  the  current  density 
increases. 

►  None  of  these  mechanisms  showed  a  downward  bent  in  the  polarization  curve. 
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Polarization  curves  for  the  porous  lithiated  NiO  cathode  are  very  often  reported  with  a  linear  slope  over 
a  wide  potential  range.  However,  the  MCFC  behaviour  at  higher  oxidant  utilization,  when  the  mass 
transfer  becomes  dominant,  is  mostly  overlooked.  Therefore,  in  this  study,  the  two  most  common 
cathode  mechanisms  are  utilized  to  compare  their  prediction  capabilities  at  higher  current  densities.  This 
is  performed  by  means  of  a  three-dimensional,  non-isothermal  mathematical  model  which  is  developed 
by  employing  volume-averaged  equations.  As  an  extension  to  previous  studies,  the  presented  model  also 
considers  the  potential  and  current  density  variation  in  both  solid  electrode  and  liquid  electrolyte  phases. 
In  essence,  this  model  is  a  set  of  partial  differential  equations  including  conservation  of  mass, 
momentum,  gaseous  species,  energy,  electronic  potential  and  ionic  potential  that  are  solved  using  a  finite 
volume  method.  In  brief,  both  peroxide  and  superoxide  mechanisms  predict  the  linear  region  of  the 
polarization  curve  accurately.  However,  none  of  these  mechanisms  showed  a  downward  bent  in  the 
polarization  curve.  A  positive  exponent  for  the  oxygen  mole  fraction  is  essential  in  obtaining  the 
downward  bent  knee  at  high  current  densities  which  is  in  contrast  to  what  has  been  reported  in  the 
literature  to-date. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Molten  carbonate  fuel  cells  (MCFCs)  are  deemed  to  be  one  of  the 
most  promising  energy  conversion  devices  that  generate  electricity 
and  heat  directly  from  oxidation  of  a  gaseous  or  gasified  fuel.  Due  to 
the  high  efficiency,  high  temperature,  quiet  operation,  flexibility  on 
fuel  and,  most  importantly,  zero  emissions,  MCFCs  have  become  an 
attractive  emerging  technology  for  stationary  co-generation  of  heat 
and  power.  Many  efforts  have  been  made  in  the  past  few  decades  to 
facilitate  the  commercialization  of  molten  carbonate  fuel  cells  and 
some  of  the  technological  concerns  are  addressed  in  [1]. 
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Nonetheless,  higher  efficiencies  and  sufficiently  elongated  oper¬ 
ating  lifetimes  with  the  least  performance  degradation  have 
received  much  attention  in  research  [1,2]. 

Generally,  adequate  models  of  porous  electrodes  would  aim  to 
enable  researchers  to  illustrate  precisely  the  relation  between  the 
MCFC  performance  and  the  structure  of  pores.  This  has  been 
a  major  subject  of  debate  among  researchers  and  several  models 
have  been  developed  to  simulate  the  behaviour  of  porous  gas- 
diffusion  electrodes.  Furthermore,  it  is  realized  that  polarization 
curves  for  the  porous  lithiated  NiO  cathode  are  very  often  reported 
with  a  linear  slope  over  a  wide  potential  range  while  unit  cell 
behaviour  at  higher  oxidant  utilization,  when  the  mass  transfer 
becomes  dominant,  is  mostly  overlooked.  Therefore,  the  goal  of  this 
paper  is  to  develop  a  three-dimensional  mathematical  model  and 
investigate  polarization  characteristics  at  higher  oxidant 
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Nomenclature 

^1,  ^2,  ^3 

anode  reaction  order 

Ti,  72 

anode  reaction  order 

Ai 

electrode  active  surface  area  (m2  m-3) 

T 

tortuosity 

A 

electrode— electrolyte  interface  are  area  (m2) 

e 

porosity 

cp 

specific  heat  (J  kg-1  K-1) 

6 

electrolyte  filling  degree 

D 

mass  diffusivity  (m2  s-1) 

V 

overpotential  (V) 

DP 

pore  diameter  (m) 

P 

dynamic  viscosity  (kg  m-1  s-1) 

£eq 

equilibrium  electric  potential  (V) 

V 

species  stoichiometric  coefficient  of  the  reaction 

A 

apparent  activation  energy  in  Eq.  (35)  (K-1) 

P 

density  (kgm-3) 

F 

Faraday’s  constant,  96,485  (Cmol-1) 

<P 

electric  potential  (V) 

*0 

exchange  current  density  (Am-2) 

a 

electric  conductivity  of  the  solid  phase  (S  m-1) 

i° 

lo 

reference  exchange  current  density  (Am-2) 

K 

electric  conductivity  of  the  electrolyte  phase  (S  m-1) 

j 

local  current  density  (Am-2) 

K0 

pre-exponential  factor  in  Eq.  (35)  (Sm-1) 

J 

average  current  density  (Am-2) 

R 

volumetric  current  density  (A  m-3) 

Subscripts  and  superscripts 

k 

thermal  conductivity  (W  m-1  K-1) 

a 

anode 

I< 

permeability  (m2) 

age 

anode  gas  channel 

M 

molecular  weight  (g  mol-1) 

c 

cathode 

n 

number  of  electrons 

ege 

cathode  gas  channel 

P 

static  pressure  (Pa) 

e 

electrolyte  phase 

Pi,  P2,  P3 

concentration  exponents  in  anode  reaction 

g 

gas  phase 

9i,  Q2,  93  concentration  exponents  in  cathode  reaction 

in 

inlet 

R 

universal  gas  constant  (8.314  J  mol-1  K-1) 

i 

ith  component 

S 

molar  entropy  (J  mol-1 1<-1) 

j 

jth  species 

S 

source  terms 

m 

mass  equation 

t 

time  (s) 

out 

outlet 

T 

temperature  (I<) 

s 

solid  phase 

gas  velocity  (ms-1) 

T 

Energy  equation 

X 

molar  fraction 

u 

momentum  equation 

Y 

mass  fraction 

(Pe 

electronic  charge  equation 

(Pc 

carbonate  ion  charge  equation 

Greek  letters 

eff 

effective 

a 

transfer  coefficient 

ref 

reference  state 

utilization.  Furthermore,  it  is  desired  to  compare  the  prediction 
capabilities  of  the  two  most  common  cathode  reaction  mecha¬ 
nisms,  namely  the  peroxide  and  superoxide  mechanisms. 

The  simplest  electrode  model,  the  so-called  Simple  (or  Flooded) 
Pore  Model,  was  introduced  by  Austin  et  al.  [3]  in  1965.  The 
problem  associated  with  this  primitive  model  was  the  very  poor 
performance  prediction  due  to  severe  mass  transfer  limitation. 
Since  then  researchers  have  tried  to  introduce  corrections  to  this 
model  which  led  to  the  Thin  Film  Model  [4,5]  and  the  Finite- 
Contact- Angle  Meniscus  Model  [6]  which  are  extensions  of  the 
simple  pore  model  and  account  for  variations  in  the  wetting 
tendency  of  the  electrolyte.  In  a  flooded  pore  with  finite-angle 
meniscus  the  current  is  predominantly  concentrated  in  a  small 
portion  of  the  pore  wall.  Whereas  in  a  film-covered  pore,  the 
electrochemical  reaction  still  has  a  propensity  to  be  concentrated  in 
the  part  of  the  film  which  is  closest  to  the  bulk  electrolyte.  Later  on, 
the  migration  of  reactants  on  the  surface  of  the  electrode  was 
incorporated  by  Iczkowski  [7]. 

Bearing  in  mind  all  these  cited  studies  a  more  realistic  porous 
electrode  model  ought  to  consider  a  spectrum  of  pore  sizes.  For  that 
reason,  several  dual-porosity  models  were  introduced  [8-10].  The 
Standard  Agglomerate  Model,  presented  by  Giner  and  Flunter  [11], 
presumes  an  idealized  electrode  in  which  the  pores  are  divided  into 
two  forms.  The  micro-pores  are  assumed  to  be  completely  flooded 
with  electrolyte,  while  the  macro-pores  are  thought  to  enclose  only 
gas.  The  model  has  been  reasonably  successful  in  predicting  the 
performance  of  fuel  cells.  The  agglomerate  model  characterizes 
electrode  structure  more  satisfactorily  than  the  thin  film  model 
[12].  Flowever,  the  anode  and  cathode  have  different  wetting 


characteristics.  Wilemski  [13]  proposed  individual  porous  elec¬ 
trode  models  for  the  anode  and  cathode  of  a  molten  carbonate  fuel 
cell.  In  this  model  all  electrochemical  activity  is  assumed  to  take 
place  on  film-covered  walls  of  the  larger  gas  filled  pores.  Smaller 
flooded  pores  were  treated  as  electrochemically  inert.  The  model 
showed  good  agreement  with  experimental  data.  Nevertheless,  it 
oversimplifies  pore  structure,  and  requires  measured  values  for 
film  areas  and  thicknesses,  forcing  these  parameters  to  be  treated 
empirically.  Afterwards,  Kunz  et  al.  [14]  modified  the  conventional 
agglomerate  model  and  developed  a  homogeneous  model  by 
calculating  the  effective  agglomerate  diameter,  porosity  and 
tortuosity  based  on  the  electrode’s  pore  spectrum  and  electrolyte 
content.  Many  researchers  have  used  this  theory  and  developed  it 
further.  Jewulski  and  Suski  [15]  proposed  an  isotropic  steady  state, 
one-dimensional  model  for  the  porous  anode  of  an  MCFC  which 
required  the  thickness  of  the  electrolyte  film  in  the  pores  as  the 
only  adjustable  parameter.  Later  on  Jewulski  [16]  applied  this 
model  for  the  porous  cathode  as  well.  Yuh  and  Selman  [12] 
developed  a  steady  state,  two-dimensional  dual-porosity 
agglomerate-type  model  for  the  porous  electrodes.  The  proposed 
model  involves  a  more  complicated  expression  instead  of  the 
Ohm’s  law  to  include  ionic  migration  in  the  melt.  They  employed 
this  model  to  predict  the  electrode  performance.  Lee  et  al.  [17] 
indicated  that  the  resulting  values  of  the  fit  parameters  depend 
strongly  on  the  choice  of  the  agglomerate  radius  (or  slab  width),  the 
film  thickness,  and  the  electrolyte  conductivity.  These  radii  and  film 
thicknesses  are  very  difficult  to  determine  meaningfully.  Further¬ 
more,  because  of  its  geometric  restrictions,  it  is  difficult  to  incor¬ 
porate  the  electrolyte  filling  degree  into  the  agglomerate  model. 
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Another  problem  of  the  agglomerate  model  is  that  the  model  does 
not  accurately  predict  the  optimal  degree  of  filling  [18]. 

However,  all  the  previous  stated  models  are  based  on 
a  continuum  approach  to  modeling.  Prins-Jansen  et  al.  [19]  chose  to 
use  a  more  fundamental  approach  named  Agglomerate-Like  Model 
which  is  based  on  averaging  of  equations  that  describe  the 
processes  on  the  micro-scale.  This  model  eliminates  the  important 
drawback  of  the  preceding  agglomerate  model  caused  by  geometric 
assumptions  and  restrictions.  Unlike  the  agglomerate  model,  the 
new  model  is  suitable  for  studying  three-dimensional  and  aniso¬ 
tropic  problems,  and  incorporates  the  degree  of  electrolyte  filling. 
Fontes  et  al.  [20]  used  a  steady  state  two-dimensional  pseudo- 
homogeneous  model  for  the  three-phase  system  of  the  cathode  and 
involved  the  polarization  curves  from  the  heterogeneous  agglom¬ 
erate  model  as  local  source  function.  They  showed  that  the  geom¬ 
etry  of  the  gas  holes  in  the  current  collector  can  influence  the 
performance  of  the  electrode  due  to  the  non-uniform  lateral 
distribution  of  the  current  production.  They  [21]  also  used  the 
homogeneous  and  heterogeneous  agglomerate  model  for  one  and 
two-dimensional  calculations,  respectively,  to  explain  the  influence 
of  gas  phase  mass  transfer  on  the  performance  of  the  porous  MCFC 
cathodes.  The  results  showed  that  for  standard  gas  composition  and 
normal  operating  current  density,  the  effect  of  gas-phase  diffusion 
is  small.  In  addition,  for  low  oxygen  partial  pressures,  the  influence 
of  mass  transfer  limitations  appears  to  be  large  even  at  low  current 
densities,  which  makes  kinetic  studies  very  difficult.  In  a  newer 
model,  called  Electrochemical-Potential  Model  [22],  the  electro¬ 
chemical  potentials  for  individual  species  are  combined  to  define 
component  potentials  which  are  separated  by  the  slow  chemical 
and/or  electrochemical  reactions.  Then,  the  reaction  rates  for  the 
slow  reactions  are  assumed  to  be  proportional  to  the  differences  in 
these  component  potentials.  Fehribach  et  al.  [23]  employed  this 
model  for  the  peroxide  mechanism  describing  the  electrochemistry 
of  an  MCFC  cathode.  Fehribach  and  Hemmes  [24]  compared  the 
polarization  losses  associated  with  the  various  diffusion-reaction- 
conduction  processes  in  MCFC  cathodes.  They  estimated  each 
type  of  polarization  loss  in  terms  of  component  electrochemical 
potentials.  The  main  advantage  of  the  component-potential 
approach  is  that  it  simplifies  both  the  analysis  and  the  computa¬ 
tions.  Also,  Subramanian  et  al.  [25,26]  employed  the  three-phase 
homogeneous  model  of  Prins-Jansen  et  al.  [19]  and  reported  the 
performance  analysis  results  based  on  a  one-dimensional  model. 

Despite  all  of  these  remarkable  efforts  in  developing  porous 
electrode  mathematical  models,  temperature  variation  effects  are 
still  overlooked.  Moreover,  hydrodynamics  of  gas  flow  in  the  gas 
channels  are  disregarded.  As  well,  the  effect  of  convective  mass  flux 
is  passed  over.  More  to  the  point,  three-dimensional  studies  based 
on  a  volume-averaging  technique  seems  to  be  absent  in  the  liter¬ 
ature.  Last  but  not  the  least,  unit  cell  behaviour  at  extreme  gas 
utilization  is  rarely  reported. 

As  a  consequence,  the  foremost  intention  of  this  study  is  to 
present  a  three-dimensional,  non-isothermal  model  by  employing 
the  volume-averaged  equations  used  by  Prins-Jansen  et  al.  [27].  As 
a  modification,  instead  of  using  a  single  equation  for  electric 
potential,  electronic  and  ionic  potential  equations  are  introduced 
and  solved  separately.  Hence,  the  presented  model  considers  the 
potential  and  current  density  variation  in  both  solid  electrode  and 
liquid  electrolyte  phases.  As  another  contribution,  gas  channels  are 
considered  explicitly  in  our  mathematical  modeling  and 
convection-diffusion  mechanisms  are  taken  into  account.  Another 
modification  is  performed  by  considering  the  degree  of  electrolyte 
filling  of  unit  cell  components. 

Fundamentally,  this  model  is  a  set  of  partial  differential  equa¬ 
tions  including  conservation  of  mass,  momentum,  gaseous  species, 
energy,  electronic  potential  and  ionic  potential.  This  mathematical 


model  will  be  implemented  for  our  future  transient,  three- 
dimensional  unit  cell  and  stack  level  simulations.  However,  in 
this  research  the  model  is  utilized  to  offer  an  improved  insight  into 
the  polarization  characteristics  of  molten  carbonate  fuel  cells.  To 
accomplish  this,  the  equations  are  discretized  and  solved  with 
ANSYS  FLUENT  [28],  a  finite  volume-based  commercial  software. 
Meanwhile,  the  programming  language  C  is  utilized  to  develop 
a  complementary  module  in  which  the  MCFC  model  details  are 
integrated.  The  developed  mathematical  model  is  then  used  to 
predict  unit  cell  behaviour  at  high  cathode  gas  utilizations. 

2.  Problem  formulation 

Fig.  1  schematically  demonstrates  the  three-dimensional  phys¬ 
ical  domain  of  an  MCFC  divided  into  five  sub-domains:  the  anode 
gas  channel  (AGC),  anode,  electrolyte,  cathode  and  cathode  gas 
channel  (CGC).  The  fuel  gas  which  is  a  gaseous  mixture  of  hydrogen, 
water  vapour  and  carbon-dioxide  enters  the  AGC  and  diffuses 
through  the  porous  anode  where  hydrogen  molecules  are  subjected 
to  the  hydrogen  oxidation  reaction  (HOR).  During  this  electro¬ 
chemical  reaction,  hydrogen  combines  with  carbonate  ions,  ending 
with  water  vapour  and  carbon-dioxide  as  follows: 

H2  +  CO^  <-►  H20  +  C02  +  2e~  ( 1 ) 

The  released  electrons  migrate  through  an  external  circuit, 
create  electricity  and  return  to  the  cell  through  the  cathode.  On  the 
cathode  side,  a  mixture  of  oxygen,  carbon-dioxide  and  nitrogen 
enters  the  cathode  gas  channel  and  diffuses  through  the  porous 
cathode,  where  the  oxygen  reduction  reaction  (ORR)  takes  place. 
Oxygen  is  reduced  to  carbonate  ions  by  combining  with  carbon- 
dioxide  and  the  electrons  coming  from  the  external  circuit  as 

lo2  +  C02  +  2e-~COf  (2) 

The  carbonate  ions  formed  at  the  cathode  move  through  the 
electrolyte  toward  the  anode,  carrying  the  electric  current  and 
completing  the  carbon-dioxide  circuit.  The  structural  parameters  of 
the  simulated  MCFC  are  summarized  in  Table  1. 

A  comprehensive  MCFC  model  needs  to  consider  the  transport  of 
the  multi-component  gas  species  in  gaseous  and  liquid  phases, 
electrochemical  and  chemical  reaction  kinetics,  heat  generation, 
heat  transfer,  transport  of  the  electrons  and  carbonate  ions,  and 
porous  electrode  effects.  These  processes  occur  in  void  volumes, 
liquid  phase,  solid  phase  and  at  triple-phase  boundaries  (TPB). 
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Table  1 

The  structural  parameters  of  the  simulated  MCFC. 


Parameter 

Value 

Cell  length,  mm 

100 

Anode  gas  channel  height,  mm 

2.0 

Anode  gas  channel  width,  mm 

2.0 

Cathode  gas  channel  height,  mm 

2.0 

Cathode  gas  channel  width,  mm 

2.0 

Anode  height,  mm 

0.7 

Anode  width,  mm 

4.0 

Cathode  height,  mm 

0.6 

Cathode  width,  mm 

4.0 

Electrolyte  height,  mm 

1.0 

Porosity  of  anode,  ea  [31] 

0.52 

Porosity  of  cathode,  ec  [31] 

0.62 

Without  losing  the  generic  physical  characteristics,  every  numerical 
simulation  is  conceived  and  developed  based  on  a  set  of  assump¬ 
tions  motivated  by  a  lack  of  experimentally  evaluated  physical 
parameters.  Likewise,  the  following  assumptions  are  made  for  this 
model:  (i)  the  chemical  species  obey  the  ideal  gas  law  and  are 
ideally  mixed;  (ii)  the  porous  anode  and  cathode  are  homogeneous; 
(iii)  the  effects  of  gravity  are  negligible;  (iv)  the  anodic  and  cathodic 
electrochemical  reactions  take  place  at  the  three-phase  boundaries 
inside  the  electrodes;  (v)  in  anode  and  cathode  sub-domain,  the  gas 
mixture  and  solid  component  are  in  a  thermal  equilibrium  state 
(have  identical  temperature)  and  hence  one  equation  with  effective 
heat  conductivity  can  describe  the  heat  transport  process;  (vi)  both 
anodic  and  cathodic  electrochemical  reactions  follow  the 
Butler— Volmer  equation;  (vii)  any  change  in  the  concentration  of 
carbonate  ions  inside  the  electrolyte  is  negligible. 

2.1.  Governing  equations 

This  section  exhibits  the  general  form  of  governing  equations 
implemented  to  model  the  MCFC.  In  contrast  to  the  approach  that 
employs  separate  differential  equations  for  different  sub-domains, 
in  this  study,  the  ‘single-domain  approach’  is  utilized.  This 
approach  considers  a  single  set  of  governing  equations  for  all  sub- 
domains.  However,  for  each  sub-domain  model  input  parameters 
(diffusivities,  conductivities,  etc.)  are  specified  separately.  As  a  result, 
no  interfacial  conditions  are  required  to  be  specified  at  internal 
boundaries  between  various  sub-domains.  Prior  to  describing  each 
specific  governing  equation,  it  is  worthwhile  to  point  out  that  each 
phenomenon  can  be  described  with  a  separate  partial  differential 
equation  which  comprises  a  transient  term,  diffusion  term, 
convection  term  and  a  source  term.  By  taking  all  these  terms  into 
account,  a  general  equation  can  be  derived  in  the  form 

^(p0)  +  V-(pT?©-reV0)  =  (3) 

where  8  is  1,  u,  Y,  ft,  4>eie  and  <\>ion  in  the  continuity,  momentum, 
species,  energy,  electronic  and  ionic  potential  equations,  respec¬ 
tively.  T©  and  5©  are  the  diffusion  coefficient  and  source  terms 
which  have  consistent  units. 

An  appropriate  approach  is  needed  to  model  the  porous  elec¬ 
trodes.  In  general,  the  agglomerate  model  divides  the  porous  media 
into  a  number  of  micro-pores  and  macro-pores,  followed  by  aver¬ 
aging  the  two-phase  equations  over  the  micro-porous  regions 
which  makes  it  difficult  to  justify  the  model  rigorously.  Essentially, 
in  this  approach  a  few  assumptions  concerning  pore  structure 
(agglomerate  radius  and  electrolyte  film  thickness)  are  unavoid¬ 
able.  Hence,  in  this  study  a  more  realistic  approach  commonly  used 
in  porous  media  problems,  the  so-called  volume  averaging,  is 
employed.  However,  according  to  porous-media  theory,  it  must  be 
feasible  to  define  representative  control  volumes  (unit  cells)  for  the 


volume-averaging  process  to  be  meaningful.  The  size  of  a  unit  cell 
must  be  chosen  with  the  intention  that  a  change  in  the  size  and/or 
position  of  the  cell  has  an  insignificant  effect  on  the  porosity  of  the 
cell.  This  implies  that  it  must  be  considerably  larger  than  the  length 
scale  of  micro-porosity,  but  much  smaller  than  the  scale  on  which 
significant  changes  in  macroscopic  quantities  arise.  In  this  model, 
all  three  phases  are  taken  into  account.  By  doing  so,  a  large  enough 
unit  cell  can  be  defined,  and  this  makes  it  possible  for  a  unit  cell  to 
fulfill  the  requirements  for  averaging.  In  contrast  to  the  agglom¬ 
erate  model,  the  new  model  is  suitable  for  studying  three- 
dimensional  and  anisotropic  problems,  and  integrating  the 
degree  of  electrolyte  filling.  In  actual  fact,  this  model  is  based  on  the 
basic  mass  and  current  balances  at  the  micro-scale  which,  subse¬ 
quently,  are  averaged  (homogenized)  across  all  three  phases  (solid/ 
gas/liquid)  of  the  electrode  to  yield  the  macro-scale  level  equations. 
In  the  following  sections,  the  macro-scale  governing  equations  are 
illustrated.  Readers  are  referred  to  Prins-Jansen  et  al.  [27]  for  more 
details  on  the  averaging  technique. 

■  Conservation  of  mass 

To  begin  with,  the  transport  of  any  gas  species  has  to  satisfy  the 
conservation  of  mass  and  momentum.  The  total  mass  gain  in  the 
anode  is  equal  to  the  mass  loss  in  the  cathode.  This  can  be  justified 
by  considering  the  production  and  consumption  of  carbon-dioxide. 
Clearly,  for  each  mole  of  CO2  produced  in  the  fuel  flow,  a  mole  of 
CO2  is  consumed  in  the  oxidant  flow  due  to  conservation  of  electric 
charges  (Eqs.  (1)  and  (2)).  Therefore,  the  conservation  of  mass  is 
written  as  [29]: 

—  (eeffPg)  +  V-  (pg T?g)  =  Sm  (4) 

where  pg  and  Ttg  are  the  superficial  values  of  the  gas  mixture 
density  and  velocity,  respectively.  Sm  (kgm-3s-1)  is  the  mass 
source  term  which  has  different  values  depending  on  the  cell  sub- 
domain.  The  gas  mixture  density  (kgm-3)  is  calculated  based  on 
the  ideal  gas  law  [29]: 

%-p«(CT  £§)  <5> 

where  Pg  is  the  superficial  value  of  the  gas  pressure  (Pa),  T  the 
temperature  (I<),  R  the  universal  gas  constant  (Jkmol-1  K-1).  In 
addition,  Yz  and  Mi  are  mass  fraction  and  molecular  weight 
(kgkmol-1)  of  species  i,  respectively. 

It  is  crucial  to  indicate  that  the  actual  volume  fractions  of  the  gas 
in  porous  anode  and  cathode  are  less  than  the  electrode  porosity 
(e).  This  is  a  result  of  the  volume  percentage  occupied  by  the 
electrolyte,  namely  the  electrolyte  filling  degree  (0),  which  is 
present  in  liquid  form.  Hence,  an  effective  porosity  (eeff)  is  defined 
and  implemented  in  the  governing  equation  as  well  as  in  consti¬ 
tutive  laws 

eeff  =  e(\  -  6)  (6) 


■  Conservation  of  momentum 

The  general  form  of  the  conservation  of  momentum  equation 
can  be  written  as  [29]: 

Of  (^egps Ug)  +  v  •  (j^ffj2ps^g^gSj  =  ~VPg  +  v  •  (*)  +  su  (7) 
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where  Su  is  the  momentum  volumetric  sink  term  (Pa  m~3)  which  is 
zero  in  gas  channels.  It  is  determined  in  the  porous  electrodes  using 
Darcy’s  equation  [29].  This  momentum  sink  contributes  to  the 
pressure  gradient  in  the  porous  electrodes,  creating  a  pressure  drop 
that  is  proportional  to  the  fluid  velocity.  For  homogeneous  porous 
media: 


S 


u 


_^Lt J 

I<e ff  U  S 


(8) 


In  this  equation  /jtg  is  the  dynamic  viscosity  of  the  ideal  gas 
mixture  (kgm-1  s-1)  and  it  is  calculated  based  on  kinetic  theory 
[29]  as 


■  Conservation  of  gaseous  species 

To  describe  the  chemical  species  transport,  the  general  form  of 
the  conservation  equation  including  both  convection  and  diffusion 
terms,  is  considered: 

(eefVi)  +  v'  (  -  PgD&VYi)  +  v'  {Pg uSYi)  =  si  04) 

where  i  represents  H2,  H20  and  C02  in  anode  and  02,  C02  and  N2  in 
cathode.  Therefore,  there  are  mainly  five  species  to  be  considered  in 
this  study.  In  order  to  facilitate  the  solution  procedure  only  the  first 
four  of  the  five  species  can  be  numerically  resolved,  leaving  the  fifth 
one  to  be  determined  through  the  following  equation: 


EWu 

j 


in  = 


[l  +  (ft/ Mj)  05 


25 


0.5 


[8  (1+Mf/M,-)] 


0.5 


(9) 


(10) 


where  i  and  j  represent  different  species.  X;  is  the  mole  fraction  of 
species  i.  Furthermore,  /<eff  is  the  effective  permeability  of  the 
porous  media  which  depends  on  the  relative  permeability,  Kr,  and 
the  intrinsic  permeability,  JQ,  through  the  following  equation  [30]: 

Kefr  =  KrI<i  (11) 

The  relative  permeability  is  determined  by  [30]: 

Kr  =  (\-d?  (12) 


EY>-  = 1  05) 

i 

Flowever,  the  mechanism  of  species  transport  in  gas  channels 
and  porous  electrodes  are  not  identical.  In  gas  channels,  no  elec¬ 
trochemical  reaction  exists  and  the  multi-component  gas  species 
transport  occurs.  Therefore,  in  Eq.  (14)  the  species  mass  source 
term,  S/,  is  zero.  is  the  effective  mass  diffusion  coefficient  of 
species  i  in  the  gas  mixture  and  is  determined  by  [29]: 


D?ff  = 


1  -Xi 


£  (Xj/Df) 


(16) 


where  D|ff,  is  the  effective  binary  mass  diffusion  coefficient  of 
species  i  in  species  j  which  is  calculated  by 


Df  =  D 


j  pref  £eff 
y’pref  p  T 


(17) 


Here,  different  values  are  used  for  the  exponent.  The  widely 
used  cubic  correlation  is  empirical  and  comes  from  sand/rock- type 
porous  media  with  a  typical  porosity  of  0.1 -0.4.  Nonetheless,  it  is 
suggested  to  be  between  4.0  and  5.0  for  porous  materials  with 
porosities  over  0.6  [30].  For  MCFCs,  it  is  recommended  to  be  1  [31  ]. 
Perceptibly,  a  combination  of  Eqs.  (11)  and  (12)  depicts  that  if  the 
local  pore  volume  of  the  anode  or  cathode  is  fully  saturated  with 
liquid  electrolyte,  the  gas  permeability  will  become  zero,  resulting 
in  an  infinite  (negative)  value  for  the  momentum  source  term.  The 
intrinsic  permeability  is  an  intensive  property.  It  is  a  measure  of  the 
ability  of  the  porous  material  to  allow  fluids  to  pass  through  it  and 
is  a  function  of  the  material  structure  only  (and  not  of  the  fluid), 
and  explicitly  distinguishes  the  value  from  that  of  relative  perme¬ 
ability.  It  is  hard  to  find  values  for  permeability  in  the  literature  for 
MCFCs.  Findlay  [32]  utilized  the  values  for  the  permeability  in 
a  proton  exchange  membrane  fuel  cell  (PEMFC)  which  is 
1.9  e  12  m2.  On  the  other  hand,  the  Carman-Kozeny  relation  for  an 
aggregated  bed  of  spheres  [33,34]  can  be  used  to  estimate  the  value 
of  the  permeability: 


I<  = 


dp 

150 


Jfl 

(1  -  £eff) 


3 


(13) 


where  Dp  is  the  pore  diameter  of  the  porous  material  which  is 
normally  between  8  and  12  pm  for  MCFC  electrodes  [32,33].  By 
substituting  this  value  in  Eq.  (13)  and  calculating  the  permeability, 
the  value  used  by  Findlay  [32]  can  be  justified.  Also,  in  a  CFD  model 
developed  by  Jiao  [35],  the  pore  diameter  of  the  catalyst  layer  is 
reported  to  be  24  pm  and  the  permeability  is  estimated  to  be 
6.2  e-12  m2.  This  obviously  shows  that  1.9  e-12  m2  is  applicable  for 
MCFCs  which  is  also  reported  by  Promislow  et  al.  [36]. 


where  Dy  is  the  bulk  binary  diffusivity  at  the  reference  temperature 
(Tref)  and  reference  pressure  (P1^).  Also,  t  is  the  tortuosity  of  the 
porous  material  which  is  frequently  estimated  by  the  following 
Bruggemann  correlation  in  fuel  cell  modeling  [30]. 

t = (,effr°'5  ns) 

It  is  apparent  that  in  gas  channels  the  porosity  is  equal  to  one. 
The  species  transport  mechanism  of  porous  electrodes  is  a  more 
complex  scenario.  Fig.  2  demonstrates  a  closer  view  of  the  porous 
electrode  morphology.  According  to  this  figure,  each  volume¬ 
averaging  cell  encloses  the  solid  electrode,  gas  mixture  and  liquid 
electrolyte.  Mass  transport  occurs  in  the  liquid  and  gas  phases. 
Precisely,  reactants  diffuse  through  the  gaseous  mixture  and  then 
transfer  to  the  molten  electrolyte  so  as  to  reach  the  triple-phase 
boundary  where  the  electrochemical  reaction  takes  place.  Hence, 
two  equations  can  be  written  for  species  transport  in  gas  and  liquid 


Fig.  2.  Schematic  of  porous  electrode  (left)  and  volume-averaging  cell  enclosing  the 
triple-phase  boundaries  (right). 
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phases  at  the  micro-scale  level,  considering  the  fact  that  electro¬ 
chemical  reactions  take  place  in  the  liquid  phase  only.  Now,  by 
defining  a  control  volume  resembling  Fig.  2,  the  two  phases  (the  gas 
and  the  electrolyte)  could  be  effectively  combined.  Clearly,  the 
physically  observable  quantities  of  interest  (concentration)  occur 
on  a  much  larger  macro-scale  where  micro-scale  equations  are  not 
applicable.  Rather,  these  micro-scale  equations  can  be  averaged 
using  theorems  from  porous-media  theory  [37].  By  doing  so,  it 
introduces  average  concentrations  defined  across  both  phases, 
represented  by  Eq.  (14).  Nevertheless,  there  is  a  major  apprehen¬ 
sion  in  regards  to  the  diffusive  flux  terms  which  is  related  to  the 
diffusion  coefficients.  It  is  said  to  be  evaluated  by  a  volume  fraction- 
based  average  over  the  gas  and  liquid  phase  diffusivities  [27].  In  this 
study  it  is  implemented  as 


mechanisms  are  summarized  in  Table  2.  The  overpotential,  rj,  is 
defined  as 

Va  =  <t>s  ~  </>e  (24) 


Vc  =  <t>s  ~  <t>e  ~  £eq  (25) 

here,  <\>s  and  <\>e  are  solid  phase  and  electrolyte  phase  potentials. 
Moreover,  Eeq  is  the  potential  difference  between  solid  and  elec¬ 
trolyte  phase  potentials  in  equilibrium,  i.e.  when  no  current  is 
generated,  and  is  defined  using  the  Nernst  equation  [39]: 


Eeq  —  Eq  + 


RT.  (PH2,aPc02,cPo25c) 
nF  {  Pc02,aPH20,a  ) 


(26) 


Deff  = 


£(1  -  6)/- 


\-Xi 


£  (*//d|) 

ij*i 


e6/D\ 


E0  =  1.2723  -  2.7645  x  10“4  T  (27) 

Now,  the  production  and  consumption  of  various  species 
involved  in  the  anode  and  cathode  reactions  can  be  expressed  as 


where  g  and  1  corresponds  to  the  gas  and  liquid  phases  respectively. 

The  source  terms  on  the  right  hand  side  of  Eqs.  (4)  and  (14)  are 
directly  associated  with  the  electrochemical  reactions  in  anode  and 
cathode.  Over  the  past  three  decades,  several  studies  have  been 
carried  out  to  find  the  anode  and  cathode  reaction  mechanisms 
that  best  represent  the  actual  electrochemical  reaction  (e.g.  [16]). 
The  available  proposed  mechanisms  are  mostly  based  on  the 
concept  of  rate-determining  step. 

The  two  most  common  reaction  mechanisms  for  anode  are 
proposed  by  Ang  and  Sammels  [38]  and  Jewulski  and  Suski  [15]. 
Boden  et  al.  [31]  incorporated  both  mechanisms  in  their  mathe¬ 
matical  modeling  and  achieved  identical  results.  In  this  study,  the 
former  mechanism  is  employed 


vi,aMiRa 

nF 


(28) 
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VicMjRc 

nF 
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where  Vi  denotes  the  stoichiometric  coefficient  of  species  i.  Also,  the 
mass  equation  source  term  reads 


Sm  =  JJS, 

i 


(30) 


■  Conservation  of  electronic  charge 


Ra  —  4V  a  X  lo  a 
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(20) 


In  contrast,  there  are  several  reaction  mechanisms  that  have 
been  proposed  for  the  cathode  electrochemical  reaction.  In  fact,  the 
prediction  capability  of  these  mechanisms  is  rarely  discussed  in 
literature.  Therefore,  the  two  most  common  reaction  mechanisms, 
namely  the  peroxide  and  superoxide  mechanisms  [12,14]  are 
employed  and  compared  in  this  paper.  In  any  event,  the  general 
form  of  the  cathode  reaction  equation  reads 


Rc 


A/,c  x  io,c 


'  xco2  \ 
XCo2,in  ) 
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exp 


'XcoA 
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exp 


(21) 


The  governing  equation  for  the  electronic  charge  transport  in 
MCFCs  can  be  derived  by  means  of  Ohm’s  law  as  follows: 

v.(-<TeffV^s)  =  S0s  (31) 

where  creff  is  the  effective  electric  conductivity  of  the  solid  material 
which  is  estimated  based  on  the  Bruggemann  correlation.  An 
exponent  of  1.0  is  used  [31]. 

<7eff  =  cr(l  —  e)1’0  (32) 

Additionally,  S ^  denotes  the  electron  generation  or  consumption  in 
the  electrodes.  It  is  determined  by  Ra  in  the  anode  and  Rc  in  the 
cathode  and  is  zero  elsewhere. 

■  Conservation  of  ionic  charge 

Similar  to  Eq.  (31),  the  ionic  charge  is  conserved  through  the 
following  equation: 

V-(xeffV0e)  =  (33) 


where  Av,  F  and  a  are  the  electrode  active  surface  area,  Faraday 
constant  and  transfer  coefficient,  respectively,  z'o  is  the  exchange 
current  density  [31],  given  by 

!0,a  =  *8, a  (%.in)A’  (*H20,in)A2  (22) 

i0,c  =  '0.c(^O2.in)T,(XCO2.in)T2  (23) 

where  i[j  is  the  standard  exchange  current  density  and  “in”  donates 
inlet.  The  values  of  transfer  coefficients  and  exponents  for  all 


Table  2 

Reaction  orders  and  species  exponents  in  the  electrochemical  reaction  rates. 


Reaction 

mechanism 

Reaction  orders 

Concentration 

exponents 

Ang  &  Sammels 

Ai  =  0.25,  X2  =  0.25, 

pi  =0.5,  p2  =  -0.5, 

(anode) 

A3  =  0.25 

P3  =  1,P4=1 

Peroxide 

71=0.375,  72  =  -1.25 

<7i  =  0  ,q2  =  -2, 

(cathode) 

Q 3  =  0.5,  (?4=-l 

Superoxide 

7i  =  0.625,  72  =  -0.75 

<7i  =  0,  q2  = -2, 

(cathode) 

<73  =  0.75,  <74  =  -0.5 
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where  /ceff  is  the  effective  electric  conductivity  in  the  liquid  phase 
which  is  approximated  based  on  the  following  correlation  for 
electrodes  [31  ] 

Keff  =  K(e0)15  (34) 

According  to  this  study,  the  effective  conductivity  of  the  pore 
electrolyte  in  the  anode  is  0.03-93  S  m-1,  depending  on  the  degree 
of  electrolyte  filling.  However,  the  following  correlation  is  used  to 
estimate  the  ionic  conductivity  [40] 

k  =  K0exp(-Ek/T)  (35) 

In  addition,  denotes  the  carbonate  ion  generation  or 
consumption  in  the  electrodes.  It  is  determined  by  Ra  in  the  anode 
and  Rc  in  the  cathode  while  it  vanishes  elsewhere. 


current  density,  respectively.  These  parameters  are  related  to  the 
potentials  through  Ohm’s  law: 

7s  =  -<rv</>s  (42) 

7  =  KV<t>e  (43) 

One  of  the  objectives  of  the  present  article,  unlike  previous 
studies,  is  to  provide  all  required  model  input  parameters  to  help 
the  reader  have  a  better  understanding  of  the  model.  Hence,  the 
electrochemical  kinetic  parameters  are  summarized  in  Table  3.  As 
well,  Table  4  exhibits  all  physical  and  thermal  properties  of  the 
materials  involved  in  the  mathematical  modeling. 

2.2.  Boundary  conditions 


■  Conservation  of  energy 

It  is  assumed  that  the  gas  mixture  and  solid  components  of  the 
fuel  cell  are  in  a  thermal  equilibrium  state  [41].  Thus,  only  one 
energy  equation  will  be  solved  for  each  cell  region.  As  such,  the 
energy  equation  applying  to  each  individual  zone  of  the  fuel  cell 
can  be  written  as  [29]: 


at  E  M  J  +  v  •  (-fceffvr)  +  V  •  (TfgepgcpT)  =  St 

\k  =  g,s,e  ) 

(36) 

The  first  term  on  the  left  hand  side  accounts  for  all  three  phases 
available  in  each  individual  sub-domain  which  incorporates  the 
volume-averaged  values  of  material  properties.  /<eff  is  the  effective 
thermal  conductivity,  determined  by 

/<eff  =  (1  -  e)ks  +  e(l  -  d)kg  +  edke  (37) 

where  ks ,  kg  and  ke  are  the  thermal  conductivity  of  solid  material, 
gas  mixture  and  liquid  electrolyte,  respectively.  Similar  to  the 
dynamic  viscosity,  kinetic  theory  is  used  to  determine  the  thermal 
conductivity  of  the  gas  mixture  as  follows: 


kg 


Xjkj 

j 


(38) 


The  heat  generation  or  consumption  is  represented  by  the 
source  term,  Sj.  Three  kinds  of  heat  sources  in  the  electrodes  are 
considered  in  the  presented  model,  namely  the  reversible  heat 
release  during  the  electrochemical  reaction,  irreversible  or  activa¬ 
tion  heat  generation  and  ohmic  heating.  The  only  heat  source  in  the 
electrolyte  is  due  to  ohmic  heating  which  is  evaluated  by 


ST  =  f  (39) 

In  the  anode,  all  three  types  of  heat  generation  mechanisms  are 
present.  Hence,  we  have 
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Likewise,  the  source  term  in  the  cathode  is 


(40) 


Based  on  the  presented  complete  set  of  governing  equations, 
there  are  12  unknowns  which  needed  to  be  solved  for 

U,  V,  W,  P,  Yh2  ,  Yh2CU  ^C02  5  Vo2  >  ^N2  Ji 

With  the  intention  of  completing  the  fuel  cell  model  formula¬ 
tion,  stating  various  boundary  conditions  at  different  positions  is 
indispensable.  The  boundary  conditions,  for  a  computational 
domain  with  a  single  pair  of  gas  flow  channels,  are  illustrated  in 
Fig.  3.  It  was  elucidated  earlier  that  the  boundary  conditions  are 
vital  only  at  the  external  surfaces  of  the  computational  domain  due 
to  the  employed  single-domain  formulation.  These  are  no-flux 
conditions  everywhere  with  the  exception  of  electrode/channel 
interfaces. 

At  the  anode  gas  channel  inlet  (/agc)  and  cathode  gas  channel 
inlet  (Jcgc),  the  total  mass  flux,  temperature  and  gas  species 
composition  of  the  entering  gas  flow  are  specified.  Moreover,  the 
fluxes  of  electric  and  ionic  charge  are  considered  to  be  zero. 
Additionally,  considering  the  very  large  aspect  ratio  of  the  gas 
channels,  the  flow  is  assumed  to  be  fully  developed  at  the  anode 
gas  channel  outlet  (Oagc)  and  cathode  gas  channel  outlet  (Ocgc).  This 
means,  none  of  the  variables  and  respective  fluxes  vary  in  the 
normal  direction.  Likewise,  the  gas  pressure  is  specified.  A  ther¬ 
mally  insulated  no  slip  boundary  condition  is  applied  to  the  anode 
and  cathode  gas  channel  walls  including,  Wjgc,  W^gc.  The  fluxes  of 
electric  and  ionic  charge  are  set  to  zero.  For  all  other  walls  of  the  gas 
channels,  W^agc,  W^agc,  w£cgc,  and  w£cgc,  the  electronic  potential 
is  specified  while  the  fluxes  corresponding  to  the  remaining  vari¬ 
ables  are  set  to  zero.  For  the  no-flux  boundaries,  W1  and  W^,  the 
zero-flux  condition  is  assumed  for  all  variables. 

2.3.  Numerical  study 

The  computational  domain  is  divided  into  a  number  of  control 
volumes  using  ANSYS  ICEM  CFD  12.0.1.  The  Finite  Volume  based 
commercial  software,  ANSYS  FLUENT  12.0.1,  is  employed  to  dis¬ 
cretize  and  solve  the  set  of  governing  equations.  These  equations 
are  strictly  coupled  through  the  source  terms  resulting  from  the 
electrochemical  reactions.  A  set  of  under-relaxation  techniques  is 
developed  to  handle  the  divergence  difficulties.  Furthermore,  the  C 
programming  language  is  used  to  develop  a  supplementary  code  in 
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Electrochemical  kinetic  parameters. 
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In  the  gas  flow  channels  no  heat  generation  occurs.  In  the  above 
equations,  Js  and  Je  are  the  magnitude  of  electronic  and  ionic 


Parameter  Value 

Standard  exchange  current  density  of  anode,  i§  a  (A  m  2)  [38]  20-220 

Standard  exchange  current  density  of  cathode,  c  (A  m-2)  [27]  0.3— 7.0 

Active  surface  area  of  anode,  Av,a  (m2  m-3)  [31  ]  2.7E5 

Active  surface  area  of  cathode,  AVyC  (m2  m-3)  [31  ]  3.0E5 
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Table  4 

Physical  and  thermal  properties  of  various  materials. 


Parameter  Value 

Thermal  conductivity  of  anode,  ka  (W  m-1  K-1)  [42,43]  78 

Thermal  conductivity  of  cathode,  kc  (W  itT1  K-1)  [42,43]  0.9 

Thermal  conductivity  of  electrolyte,  ke  (W  itt1  K-1)  [43]  2.0 

Specific  heat  of  anode,  cp,a  (J  kg-1 I<-1)  [46]  444 

Specific  heat  of  cathode,  cpx  (J  kg-1  K_1)  [46]  44,352 

Specific  heat  of  electrolyte,  cpx  (J  kg-1  K-1)  [46]  4000 

Density  of  anode,  pa  (kg  itt3)  [46]  8220 

Density  of  cathode,  pc  (kg  m-3)  [46]  6794 

Density  of  electrolyte,  pe  (kg  m  3)  [46]  1914 

Electric  conductivity  of  anode,  ua  (S  m-1)  [20,44]  100 

Electric  conductivity  of  cathode,  ua  (S  m-1)  [20,44]  100 

Free  electrolyte  conductivity:  pre-exponential  factor,  3637 

/c0  (S  m-1)  [13] 

Free  electrolyte  conductivity:  apparent  activation  energy,  3016 

£fc(K  ’)  [13] 

Hydrogen  diffusivity  in  carbon-dioxide,  DH2-co2  (m2  s_1)  [46,47]  5.5E-5 

Hydrogen  diffusivity  in  water  vapour,  £*h2-h2o  (rn2s_1)  [46,47]  9.15E-5 

Oxygen  diffusivity  in  carbon-dioxide,  D02-co2  (m2  s  1)  [46,47]  1.4E-5 

Oxygen  diffusivity  in  nitrogen,  Dq2-n2  (m2  s_1)  [46,47]  1.8E-5 

Carbon-dioxide  diffusivity  in  water  vapour,  1.62E-5 

£>co2— h2o  (m2  s_1)  [46,47] 

Carbon-dioxide  diffusivity  in  nitrogen,  DC02-n2  (m2  s_1)  [46,47]  1.6E-5 

Hydrogen  diffusivity  in  liquid  electrolyte,  (m2  s_1)  [25]  IE-7 

Oxygen  diffusivity  in  liquid  electrolyte,  (m2  s_1)  [25,29]  3E-7 

Carbon-dioxide  diffusivity  in  liquid  electrolyte,  IE-7 

(m2s-1)  [25,29] 

Water  vapour  diffusivity  in  liquid  electrolyte,  DlH  0  (m2  s_1)  IE-7 

Nitrogen  diffusivity  in  liquid  electrolyte,  (m2  s_1)  IE-7 

Standard  entropy  change  of  anode,  ASa  (J  moF1 1<-1)  [45]  54.56 

Standard  entropy  change  of  cathode,  ASC  (J  mol-1  K_1)  [45]  -216.2 

Standard  entropy  change  of  generation  reaction,  -161.64 

AStot  (J  mol-1  K_1)  [45] 


order  to  add  several  capabilities  (e.g.  non-standard  transport 
equations,  source  terms,  temperature  dependent  properties,  etc.)  to 
ANSYS  FLUENT  12.0.1.  The  SIMPLE  algorithm  is  implemented  for  the 
coupling  between  the  pressure  and  velocity  field.  An  algebraic 
multi-grid  (AMG)  method  with  a  Gauss— Seidel  type  smoother  is 
used  to  accelerate  the  convergence.  A  strict  convergence  criterion 
with  a  residual  of  10  10  is  used  for  all  variables. 

3.  Results  and  discussion 


polarization  curve.  In  the  first  step,  the  model  is  built  based  on  the 
reported  geometry  and  operating  conditions  in  [49].  In  brief,  the 
authors  conducted  several  experiments  to  examine  the  perfor¬ 
mance  of  a  single  MCFC  with  a  planar  area  of  100  cm2.  The 
elcotrode-electrolyte  assembly  consisted  of  a  0.77  mm  Ni-Cr  alloy 
anode  and  a  0.72  mm  in  situ  oxidized  NiO  cathode.  A  62  mol  % 
Li2C03,  38  mol  %  K2CO3  eutectic  carbonate  was  utilized  as  electro¬ 
lyte,  and  LiA102  with  Al203  fiber  was  fabricated  as  matrix  material. 
At  the  anode  side,  a  gaseous  mixture  of  H2,  H20  and  C02  with  molar 
fractions  of  0.69,  0.14  and  0.17  was  used  while  in  the  cathode  side, 
a  mixture  of  02,  C02  and  N2  with  molar  fractions  of  0.15,  0.30  and 
0.55  was  fed  into  the  system.  Base  on  these  conditions  the  polari¬ 
zation  curve  at  various  operating  condition  was  obtained  [49]. 

In  this  study,  the  polarization  curve  at  atmospheric  pressure  is 
initially  used  to  verify  and  adjust  the  model  input  parameters.  The 
two  most  common  mechanisms  for  the  cathode  electrochemical 
reaction  rate,  namely  the  peroxide  and  superoxide  mechanisms,  are 
employed  to  investigate  the  validity  of  the  model.  Details  of  these 
mechanisms  will  be  discussed  in  the  following  section.  In  the  next 
step,  the  verified  model  is  put  in  to  practice  to  predict  the  polari¬ 
zation  curve  for  another  MCFC,  reported  in  [48].  This  step  serves  as 
the  validation  process.  In  conclusion,  the  verification  and  validation 
outcome  is  summarized  in  Figs.  4  and  5.  The  former  illustrates  the 
measured  [49]  and  verified  values  based  on  the  model  fine-tuning. 
In  contrast,  the  latter  demonstrates  the  accuracy  of  the  verified 
model  in  comparison  with  another  experimental  study  [48]  and 
also  a  numerical  study.  It  can  be  easily  observed  that  the  model 
prediction  have  matched  the  experimental  data  adequately.  Based 
on  these  two  strict  steps  of  verification  and  validation  processes, 
the  developed  mathematical  model  is  deemed  functional  in  the 
sense  that  the  model  addresses  the  right  problem  and  presents 
accurate  results.  Thereafter,  the  validated  model  is  engaged  in 
several  simulation  cases  to  examine  polarization  characteristics  of 
a  unit  cell  at  relatively  high  utilization  for  cathode  gas. 

3.2.  MCFC  behaviour  at  high  cathode  gas  utilization 

In  general,  the  first  region  of  a  fuel  cell  polarization  curve,  the 
so-called  activation  polarization  region,  occurs  due  to  the  slug¬ 
gishness  of  the  electrochemical  reaction  at  low  current  densities  (or 


3.1.  Model  verification  and  validation 

It  is  generally  acknowledged  that  model  validation  is  the  most 
significant  step  in  the  model  building  sequence  of  numerical 
simulations.  It  is  also  one  of  the  most  overlooked.  In  this  study,  with 
the  purpose  of  achieving  highly  accurate  results,  the  experimental 
study  by  Brouwer  et  al.  [48]  and  Lee  et  al.  [49]  are  selected  to 
investigate  the  model  trustworthiness.  In  fact,  the  mathematical 
model  is  validated  by  means  of  the  most  popular  criterion,  the 
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Fig.  3.  Boundary  conditions  for  the  MCFC  model:  (a)  front  view,  (b)  side  view. 


Fig.  4.  Model  parameter  verification:  based  on  peroxide  and  superoxide  mechanism  for 
cathode  electrochemical  reaction  rate  using  the  experimental  study  by  Lee  et  al.  [49]. 
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Fig.  5.  Model  validation:  based  on  superoxide  mechanism  for  cathode  electrochemical 
reaction  rate  using  the  experimental  and  numerical  study  by  Brouwer  et  al.  [48]. 


voltage  losses)  across  the  cell.  Boosting  the  supplied  gas  tempera¬ 
ture  provides  some  of  the  activation  energy  required  to  drive  the 
electrochemical  reaction.  However,  in  molten  carbonate  fuel  cells 
this  region  vanishes  since  the  operating  temperature  is  very  high. 
This  could  be  observed  both  in  Figs.  4  and  5.  Therefore,  the  polar¬ 
ization  curve  straight  away  starts  with  a  linear  behaviour  which  is 
due  to  the  ohmic  polarization.  This  is  the  practical  and  useful  region 
of  any  fuel  cell  which  is  normally  followed  by  a  bent  downwards. 
Conversely,  opposite  behaviour  is  observed  in  the  simulated  MCFC, 
shown  in  Fig.  6.  In  fact,  Fig.  6  implies  that  using  the  prevailing 
reaction  mechanisms,  the  linear  behaviour  could  be  achieved 
flawlessly  but  only  for  low  and  moderate  current  densities  or 
cathode  gas  utilizations.  It  is  worthwhile  to  mention  that  almost  all 
previous  experimental  and  numerical  studies  for  MCFCs  have  only 
shown  completely  linear  polarization  curves.  In  Fig.  6,  first  the 
experimental  data  [49],  which  were  used  for  the  verification 
process,  are  extrapolated  linearly  in  order  to  estimate  the  MCFC 
polarization  characteristics  at  higher  current  densities.  Then  the 
model  predictions,  based  on  the  two  most  common  cathode 
mechanisms,  are  included.  It  is  observed  that  model  predictions 
demonstrate  an  upwards  bent.  This  has  been  rarely  discussed  in  the 
literature. 

This  paper  reveals  why  the  downwards  curve  is  not  reported  in 
numerical  studies.  Noticeably,  the  issue  was  found  to  be  in  the 
porous  cathode.  In  actual  fact,  it  is  expected  that  as  the  gaseous 
reactants  are  being  consumed  to  produce  current,  the  mass  transfer 
process  becomes  dominant  in  controlling  the  fuel  cell  performance. 
In  other  words,  the  availability  of  gas  mixture  components  turns 
out  to  be  crucial  which  should  lead  to  a  drop  off  in  the  amount  of 
current  generation.  Nevertheless,  Figs.  7  and  8  show  different 
trends.  These  graphs  demonstrate  the  variation  of  local  cathode 
volumetric  current  density  (or  carbonate  ion  generation  rate)  with 
carbon-dioxide  utilization  for  various  overpotentials  using 
peroxide  and  superoxide  mechanisms.  Both  mechanisms  show  an 
exponentially  rising  tendency  (for  all  overpotential  values)  as  the 
carbon-dioxide  mole  fraction  is  decreased.  The  peroxide  mecha¬ 
nism  (Fig.  7)  indicates  the  fact  that  the  volumetric  generation  of 
current  is  increasing  up  to  a  point  where  there  is  almost  no  carbon- 
dioxide  left  in  the  cathode.  Higher  overpotentials  exacerbate  the 
situation.  The  same  growing  trend  may  be  observed  from  Fig.  8  for 


Fig.  6.  Evaluation  of  polarization  curve  linearity  between  the  model  predictions 
(superoxide  and  peroxide  mechanisms)  and  a  linear  extrapolation  of  experimental 
data  [49]. 

the  superoxide  mechanism  which  is  less  utilized  compared  with 
the  peroxide  mechanism.  However,  for  lower  values  of  over¬ 
potentials,  the  rate  of  reaction  decreases  as  the  utilization  reaches 
almost  80%.  Therefore,  because  of  the  tendency  in  reaction  rate  (or 
volumetric  current  density)  to  increase  with  falling  CO2,  it  is 
obvious  that  these  mechanisms  will  never  predict  a  downwards 
bent  in  the  polarization  curve.  This  could  be  justified  by  looking 
back  at  Eq.  (21)  which  describes  the  cathode  electrochemical 
reaction  rate  or  volumetric  current  density.  It  suggests  that  the 
negative  exponent  of  the  oxygen  mole  fraction  causes  larger  values 
at  lower  mole  fractions.  According  to  Table  2,  the  magnitude  of  this 
negative  value  is  greater  for  the  peroxide  mechanism  which  creates 
a  bigger  divergence  from  the  linear  curve  in  Fig.  6. 

In  conclusion,  the  literature  seems  to  lack  experimental  and 
numerical  data  for  polarization  curve  at  extremely  high  cathode  gas 
utilizations  (high  average  current  densities).  Investigations  have 


Fig.  7.  Variation  of  cathode  volumetric  current  density  with  carbon-dioxide  utilization 
in  various  overpotentials  using  peroxide  mechanism  for  cathode  reaction  rate. 
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Fig.  8.  Variation  of  cathode  volumetric  current  density  with  carbon-dioxide  utilization 
in  various  overpotentials  using  superoxide  mechanism  for  cathode  reaction  rate. 


always  been  conducted  at  low  and  moderate  current  densities, 
namely  the  practical  operating  region.  This  research  shows  that 
none  of  the  available  cathode  reaction  mechanisms  can  predict  the 
MCFC  behaviour  at  high  cathode  gas  utilizations.  This  issue  is 
related  to  the  current  correlations  for  cathode  reaction  rate  [12,14] 
which  characterize  only  apparent  values  without  considering  the 
true  kinetics.  Perhaps,  one  way  to  obtain  a  good  fit  to  the  experi¬ 
mental  data  at  high  current  densities  is  to  assign  a  positive  expo¬ 
nent  to  CO2  in  Eq.  (21).  Changing  this  exponent  at  very  low 
concentration  could  be  another  alternative.  In  any  case,  conducting 
further  experimental  analysis  would  be  helpful  to  check  whether 
the  prediction  in  Fig.  9  is  accurate.  In  contrast  to  the  MCFCs,  the 
downwards  bent  exists  in  the  polarization  curve  of  proton 
exchange  membrane  fuel  cells  (e.g.,  [50,51])  and  solid  oxide  fuel 
cells  ([52,53]).  They  both  have  positive  exponents  for  reaction 
orders. 
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Fig.  9.  Polarization  curve  obtained  by  using  superoxide  mechanism  for  cathode 
reaction  rate. 
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Fig.  10.  Local  current  density  (Am  2)  distribution  in  anode  and  cathode  outlet-anode 
mechanism:  Ang  and  Sammels,  cathode  mechanism:  Superoxide  (operating  voltage: 
0.28  V). 

Fig.  9  is  obtained  by  conducting  several  simulations  for  different 
cell  voltages.  Each  simulation  results  a  point  (voltage-current 
density)  in  this  figure.  By  connecting  these  points  the  polarization 
curve  is  achieved.  However,  the  code  did  not  converge  for  the 
simulations  in  which  the  cell  voltage  was  below  0.28  V. 

In  general,  at  the  end  of  the  simulation,  once  the  electrolyte 
phase  and  solid  phase  potentials  are  determined  in  the  anode, 
cathode  and  electrolyte,  the  local  current  density  is  calculated  using 
Eqs.  (42)  and  (43).  Thereafter,  the  average  current  density  of  the 
electrode  is  determined  by 


where  A  is  the  electrode— electrolyte  interface  area,  which  should 
match  for  both  anode  and  cathode.  In  the  case  that  the  simulation 
reaches  convergence,  the  calculated  average  current  density  in 
anode  and  cathode  are  identical.  As  well,  the  local  distribution  of 
current  density  in  anode  and  cathode  are  similar.  In  contrast, 
according  to  Fig.  10,  non-uniformity  in  local  volumetric  current 
density  started  to  occur  after  V  =  0.28  V  and  different  values  for 
average  current  density  were  obtained  (Fig.  11).  Apparently,  Fig.  10 
shows  that  by  lowering  the  cell  voltage  right  after  0.28  V,  which 
corresponds  to  0.79  V  drop  in  cell  voltage,  the  local  current  density 
at  the  cathode  corners  (shown  in  dotted-box)  diminishes.  In  other 
words,  it  declines  by  a  factor  of  2  comparing  to  the  anode  local 
current  density  at  the  same  spot.  This  could  be  explained  by  taking 
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Fig.  11.  Local  volumetric  current  density  (Am-3)  distribution  in  anode  and  cathode 
inlet  and  outlet-anode  mechanism:  Ang  and  Sammels,  cathode  mechanism:  Super¬ 
oxide  (operating  voltage:  0.28  V). 
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Fig.  12.  Superoxide  mechanism  prediction  for  carbon-dioxide  mole  fraction  at  cathode 
and  cathode  gas  channel  (operating  voltage:  0.28  V). 


Fig.  11  into  consideration.  It  simply  demonstrates  the  sluggishness 
of  the  reaction  rate  or  volumetric  current  density  at  the  corners.  In 
fact,  in  this  extreme  condition,  the  cathode  electrochemical  reac¬ 
tion  occurs  mainly  under  the  gas  channel  and  partially  in  a  small 
spot  close  to  the  gas  channel.  Since  the  rate  of  reaction  is  decreased, 
the  amount  of  generated  current  is  reduced.  Consequently,  non¬ 
uniformity  of  the  local  current  density  deteriorates  the  MCFC 
performance.  In  realistic  operating  condition,  the  same  current  has 
to  transfer  through  electrode-electrolyte  assembly.  However,  the 
local  current  density  is  not  necessarily  uniform.  The  surface  area 
through  which  the  current  is  transferred,  determines  the  current 
density  and  its  area-weighted  surface  integral  identifies  the  cell 
total  current  density.  Our  mathematical  model  perfectly  predicts 
this  equality  hence  this  value  was  used  to  generate  the  polarization 
curves  up  to  V  =  0.28  V.  At  this  voltage,  the  anode  and  cathode 
average  current  density  calculated  as  6343  and  5740  A  m-2  which 
means  they  are  not  identical  to  any  further  extent.  This  is  all  due  to 
the  extremely  small  molar  fraction  of  carbon-dioxide  (and  conse¬ 
quently  oxygen)  in  the  cathode  corners  shown  in  Fig.  12.  This  figure 
shows  the  local  molar  fraction  of  carbon-dioxide  in  cathode  and 
cathode  gas  channel  inlet  and  outlet.  It  is  obvious  that  carbon- 
dioxide  utilization  reaches  its  highest  value  at  the  cathode 
corners  and  the  electrochemical  reaction  rate  finds  smaller  values 
at  these  spots.  Therefore,  the  lack  of  reactants  in  high  current 
densities  (or  high  utilizations)  is  the  source  of  all  non-uniformity 
and  performance  deterioration. 

4.  Conclusions 

A  detailed  three-dimensional,  mathematical  model  of  an  MCFC 
was  presented  by  employing  volume-averaged  equations.  Two 
different  mechanisms  were  incorporated  for  the  cathode  electro¬ 
chemical  reaction  rate.  ANSYS  FLUENT  12.0.1,  a  finite  volume  based 
commercial  software  package,  was  utilized  to  solve  the  system  of 
partial  differential  equations.  The  following  concluding  remarks  are 
drawn  from  this  study: 

■  For  a  specific  case,  available  experimental  data  were  used  to 
adjust  model  input  parameter. 

■  The  verified  model  was  employed  to  check  the  validity  of  the 
results  by  conducting  a  comparison  with  another  experimental 
study  and  good  agreement  was  achieved. 

■  Both  peroxide  and  superoxide  mechanisms  predicted  the 
linear  region  of  the  polarization  curve  accurately.  This  occurs 
for  low  and  moderate  voltage  losses,  or  low  to  moderate 
cathode  gas  utilizations. 

■  Both  mechanisms  showed  a  concave  tendency  for  the  cathode 
reaction  rate  as  the  carbon-dioxide  mole  fraction  is  decreased 
when  the  current  density  increases. 


■  None  of  these  mechanisms  showed  a  downward  bent  in  the 
polarization  curve. 

■  The  negative  exponent  of  the  oxygen  mole  fraction  was  iden¬ 
tified  to  be  the  causes  of  larger  reaction  rates  at  lower  O2  mole 
fractions. 

■  Larger  negative  exponent  for  the  cathode  reaction  rate  of  the 
peroxide  mechanism  (as  compared  to  the  superoxide  mecha¬ 
nism)  created  a  larger  divergence  from  the  linearity  in  the 
polarization  curve. 

■  Using  positive  exponents  for  the  oxygen  in  the  cathode  reac¬ 
tion  rate  would  probably  result  in  obtaining  a  good  fit  to  the 
experimental  data  at  high  current  densities. 

■  At  extreme  conditions  (high  voltage  losses),  the  local  current 
density  at  the  cathode  corners  declines  by  a  factor  of  2. 

■  At  high  voltage  drops,  the  cathode  electrochemical  reaction 
occurs  mainly  under  the  gas  channel  and  partially  in  a  small 
location  close  to  the  gas  channel. 

■  The  lack  of  reactants  at  high  current  densities  (or  high  utili¬ 
zations)  is  the  source  of  all  non-uniformity  and  performance 
deterioration. 
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